Visualize Dual-Doppler Output from the Southern Great Plains Site

Imports

import xarray as xr
from distributed import Client
import hvplot.xarray
import holoviews as hv
import glob
from pathlib import Path
import numpy as np
import pandas as pd
import panel as pn
from bokeh.models.formatters import DatetimeTickFormatter
import pyart
hv.extension("bokeh")
hv.output(widget_location='bottom')
import warnings
warnings.filterwarnings("ignore")

client = Client()
client

2022-09-29 17:06:37,160 - distributed.diskutils - INFO - Found stale lock file and directory '/var/folders/bw/c9j8z20x45s2y20vv6528qjc0000gq/T/dask-worker-space/worker-k820hxo9', purging
2022-09-29 17:06:37,161 - distributed.diskutils - INFO - Found stale lock file and directory '/var/folders/bw/c9j8z20x45s2y20vv6528qjc0000gq/T/dask-worker-space/worker-n1mful2t', purging
2022-09-29 17:06:37,161 - distributed.diskutils - INFO - Found stale lock file and directory '/var/folders/bw/c9j8z20x45s2y20vv6528qjc0000gq/T/dask-worker-space/worker-swrpukpz', purging
2022-09-29 17:06:37,161 - distributed.diskutils - INFO - Found stale lock file and directory '/var/folders/bw/c9j8z20x45s2y20vv6528qjc0000gq/T/dask-worker-space/worker-a9h_h0nl', purging
2022-09-29 17:06:37,161 - distributed.diskutils - INFO - Found stale lock file and directory '/var/folders/bw/c9j8z20x45s2y20vv6528qjc0000gq/T/dask-worker-space/worker-d290v_b1', purging
2022-09-29 17:06:37,161 - distributed.diskutils - INFO - Found stale lock file and directory '/var/folders/bw/c9j8z20x45s2y20vv6528qjc0000gq/T/dask-worker-space/worker-w617ai1b', purging
2022-09-29 17:06:37,161 - distributed.diskutils - INFO - Found stale lock file and directory '/var/folders/bw/c9j8z20x45s2y20vv6528qjc0000gq/T/dask-worker-space/worker-9y_pzy9u', purging
2022-09-29 17:06:37,162 - distributed.diskutils - INFO - Found stale lock file and directory '/var/folders/bw/c9j8z20x45s2y20vv6528qjc0000gq/T/dask-worker-space/worker-ivx9z48z', purging
2022-09-29 17:06:37,162 - distributed.diskutils - INFO - Found stale lock file and directory '/var/folders/bw/c9j8z20x45s2y20vv6528qjc0000gq/T/dask-worker-space/worker-2llesuyw', purging

Client

Client-fc8f63bc-4042-11ed-a77b-520a01803a92

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:8787/status

Cluster Info

Read the Data and Break into Different Cases

All of our data has been pre-transferred from the Argonne laboratory computing facility, using the ../../data/get_data_from_lcrc.sh script.

all_cases = sorted(glob.glob("../../data/sgp_dda/*.nc"))

Determine the unique days from the files

days = []
for file in all_cases:
    days.append(Path(file).stem.split(".")[0])
unique_days = sorted(np.unique(np.array(days)))
unique_days

['20180529',
 '20180530',
 '20180531',
 '20180630',
 '20180701',
 '20180714',
 '20180715',
 '20180726',
 '20180816',
 '20180819']

Add latitude and longitude data to the grid

By default, the grid is missing latitude and longitude values. We add those here using a function from Py-ART, assuming a azimuthal equadistant projection. We only need to do this for single file, since the grid is uniform across the campaign.

ds = xr.open_dataset("../../data/sgp_dda/20180529.000800.nc")
x, y = np.meshgrid(ds.x.values, ds.y.values)
lon, lat = pyart.core.cartesian_to_geographic_aeqd(x, y, ds.origin_longitude.values[0], ds.origin_latitude.values[0])

Plot the data

We setup a loop here, looking at two main plots:

  • A timeseries of the peak vertical velocity in the midlevels (3-8 km) across

    • All vertical levels in that range

    • At each vertical level in that range

  • A domain plot (across latitude and longitude) with vertical velocity and reflectivity included

    • Note - we can use the slider to look at different vertical levels, and times

We are using hvPlot here, which enables us to produce interactive plots!

domain_plots = []
w_plots = []
for day in unique_days:
    plots = pn.Column()
    try:
        ds = xr.open_mfdataset(f"../../data/sgp_dda/{day}*")
        ds["lon"] = (("y", "x"), lon)
        ds["lat"] = (("y", "x"), lat)
        ds = ds.set_coords(["lat", "lon"])
    except RuntimeError:
        continue
    
    w_95 = np.round(ds.where(ds.w > 0).sel(z=slice(3_000, 8_000)).w.compute().quantile(.95, skipna=True).values, 2)
    
    # Adjust the long name (used in plotting) to be a bit shorter
    ds.z.attrs["long_name"] = 'height_above_ground'
    
    # Plot up our vertical velocity and reflectivity data using hvPlot, which creates an interactive plotting interface
    w_plot = ds.w.hvplot.contourf(cmap='Reds',
                                  title=f'95th Percentile Vertical Velocity {w_95} m/s',
                                  clim=(w_95, 15),
                                  x='lon',
                                  y='lat',
                                  levels=10,
                                  height=600,
                                  clabel='Vertical Velocity (m/s)',
                                 )
    z_plot = ds.corrected_reflectivity.hvplot.quadmesh(cmap='pyart_HomeyerRainbow',
                                                       clim=(-20, 60),
                                                       clabel='Reflectivity (dBZ)',
                                                       x='lon',
                                                       y='lat',
                                                       height=600,
                                                      )
    
    formatter = DatetimeTickFormatter(hours="%d %b %Y \n %H:%M UTC")
    all_levels = ds.w.sel(z=slice(3000, 8000)).max(dim=["z", "y", "x"]).hvplot.line(x='time',
                                                                                    ylabel='Peak Midlevel (3-8 km) \n Vertical Velocity in Domain (m/s)',
                                                                                    xformatter=formatter)
    each_level = ds.w.sel(z=slice(3000, 8000)).max(dim=["y", "x"]).hvplot.quadmesh(x='time',
                                                                                   ylabel='Peak Midlevel (3-8 km) \n Vertical Velocity in Domain (m/s)',
                                                                                   cmap='Reds',
                                                                                   clabel='Vertical Velocity (m/s)',
                                                                                   xformatter=formatter)
    
    w_profile_plot = (all_levels + each_level).cols(1)
    
    domain_plot = (w_plot + z_plot).cols(2)
    
    # Determine the start and end time for the title
    #time = pd.to_datetime(ds.time)
    #start_time = time[0].strftime("%b %d %Y %H:%M:%S UTC")
    #end_time = time[-1].strftime("%b %d %Y %H:%M:%S UTC")
    
    # Add a title to the subsection
    #plots+=pn.pane.Markdown(f"## {start_time} - {end_time}")
    
    # We want to add these plots to a list of all the days
    #plots+=plot
    domain_plots.append(domain_plot)
    w_plots.append(w_profile_plot)

Exception ignored in: <bound method GCDiagnosis._gc_callback of <distributed.utils_perf.GCDiagnosis object at 0x127001ea0>>
Traceback (most recent call last):
  File "/Users/mgrover/miniforge3/envs/pyart-docs/lib/python3.10/site-packages/distributed/utils_perf.py", line 189, in _gc_callback
    self._fractional_timer.start_timing()
  File "/Users/mgrover/miniforge3/envs/pyart-docs/lib/python3.10/site-packages/distributed/utils_perf.py", line 116, in start_timing
    assert self._cur_start is None
AssertionError: 
/Users/mgrover/miniforge3/envs/pyart-docs/lib/python3.10/site-packages/dask/array/reductions.py:615: RuntimeWarning: All-NaN slice encountered
  return np.nanmax(x_chunk, axis=axis, keepdims=keepdims)
/Users/mgrover/miniforge3/envs/pyart-docs/lib/python3.10/site-packages/dask/array/reductions.py:615: RuntimeWarning: All-NaN slice encountered
  return np.nanmax(x_chunk, axis=axis, keepdims=keepdims)
/Users/mgrover/miniforge3/envs/pyart-docs/lib/python3.10/site-packages/dask/array/reductions.py:615: RuntimeWarning: All-NaN slice encountered
  return np.nanmax(x_chunk, axis=axis, keepdims=keepdims)
/Users/mgrover/miniforge3/envs/pyart-docs/lib/python3.10/site-packages/dask/array/reductions.py:615: RuntimeWarning: All-NaN slice encountered
  return np.nanmax(x_chunk, axis=axis, keepdims=keepdims)
/Users/mgrover/miniforge3/envs/pyart-docs/lib/python3.10/site-packages/dask/array/reductions.py:615: RuntimeWarning: All-NaN slice encountered
  return np.nanmax(x_chunk, axis=axis, keepdims=keepdims)

View our Timeseries and Domain Plots

Now, let’s visualize this! We need to use the following to actually view the interactive plots we produced.

for plot_number in range(len(domain_plots)):
    display(w_plots[plot_number])
    display(domain_plots[plot_number])

/Users/mgrover/miniforge3/envs/pyart-docs/lib/python3.10/site-packages/dask/array/reductions.py:586: RuntimeWarning: All-NaN slice encountered
  return np.nanmin(x_chunk, axis=axis, keepdims=keepdims)
/Users/mgrover/miniforge3/envs/pyart-docs/lib/python3.10/site-packages/dask/array/reductions.py:586: RuntimeWarning: All-NaN slice encountered
  return np.nanmin(x_chunk, axis=axis, keepdims=keepdims)
/Users/mgrover/miniforge3/envs/pyart-docs/lib/python3.10/site-packages/dask/array/reductions.py:586: RuntimeWarning: All-NaN slice encountered
  return np.nanmin(x_chunk, axis=axis, keepdims=keepdims)
/Users/mgrover/miniforge3/envs/pyart-docs/lib/python3.10/site-packages/dask/array/reductions.py:586: RuntimeWarning: All-NaN slice encountered
  return np.nanmin(x_chunk, axis=axis, keepdims=keepdims)
/Users/mgrover/miniforge3/envs/pyart-docs/lib/python3.10/site-packages/dask/array/reductions.py:586: RuntimeWarning: All-NaN slice encountered
  return np.nanmin(x_chunk, axis=axis, keepdims=keepdims)